Droplet motion driven by humidity gradients during evaporation and condensation

Abstract The motion of droplets on solid surfaces in response to an external gradient is a fundamental problem with a broad range of applications, including water harvesting, heat exchange, mixing and printing. Here we study the motion of droplets driven by a humidity gradient, i.e. a variation in concentration of their own vapour in the surrounding gas phase. Using lattice-Boltzmann simulations of a diffuse-interface hydrodynamic model to account for the liquid and gas phases, we demonstrate that the droplet migrates towards the region of higher vapour concentration. This effect holds in situations where the ambient gradient drives either the evaporation or the condensation of the droplet, or both simultaneously. We identify two main mechanisms responsible for the observed motion: a difference in surface wettability, which we measure in terms of the Young stress, and a variation in surface tension, which drives a Marangoni flow. Our results are relevant in advancing our knowledge of the interplay between gas and liquid phases out of thermodynamic equilibrium, as well as for applications involving the control of droplet motion. Graphic abstract Supplementary Information The online version contains supplementary material available at 10.1140/epje/s10189-024-00426-7.

In order to numerically solve the equations of state under this framework, we make use of the Lattice Boltzmann method, a solver which considers both macroscopic and microscopic interactions.The Boltzmann equation, which arises from the Kinetic Theory of Gases [1], is a statistical description of the evolution of a fluid.The discretised form of the Boltzmann equation is what is referred to a the Lattice Boltzmann equation [2,3] To perform the Lattice Boltzmann method, a distribution function is defined, f q (x), that corresponds to the mean density of particles found at an instance in time t in a discrete position x and a velocity c q [4].
The Lattice-Boltzmann Method consists of two main steps [2,4,5], the first one is called the streaming step which describes the mean density of particles in a given position and velocity at a certain time and is given by the distribution function The second step is called the collision step, which along post-collision terms, describes how the distribution of particles has changed due to interactions in the streaming step.The distribution function associated with the collision step is expressed as follows where Ξ[f ] q is the collision operator and S q is the sources term.
During the streaming step, the particle populations interact with neighboring lattice nodes by means of the velocity c q over a unitary time step such that x → x+c q .This defines a set of velocities associated with neighboring lattice points (depending on the dimension of the domain) in the form of c q ∈ {c q } Q−1 q=0 .The result of this interaction are a Q, D-dimensional integer vectors that describe the connectivity between the lattice points.This is synthesized by the Q-D notation in the Lattice Boltzmann model.For the case presented in this work, we use a D2Q9 in the Lattice Boltzmann notation, which implies a two-dimensional system with a connectivity of 3x3 neighbourhood.equilibrium distribution function, f e q , is defined.This function defines different relaxation rates for each component of the distribution function.The collision term bears a significant amount of importance in physical terms since it enables the user to tune transport coefficients such as the viscosity [4].
The source term in Equation 2, on the other hand, can take many forms which depend on the type of sources the user wants to be considered in the simulation, for instance, defining sources of mass, forces, or stresses.
From the Navier-Stokes equation, it is possible to obtain macroscopic parameters such as the density and the velocity.By using the Lattice Boltzmann Method to numerically integrate the Navier-Stokes equations, one can define the expression of macroscopic variables in terms of the distribution function [4].The mass density flux density yields and the momentum of the density, from where one can obtain the velocity field, is defined as [6] ρu The Cahn-Hilliard equation allows us to obtain the value of the phase field.In order to compute it, a second distribution function, g q , of similar form as Equations 1 and 2 is used where the moments and collision parameters differ from f q .The zeroth moment of this distribution function defines the phase field as [6] To calculate the values of the pressure tensor field and the chemical potential, it is necessary to obtain the gradient and the Laplacian of the phase field by a 3x3 finite differences stencil.These can be approximated by using the finite difference stencils, which for the gradient yields and for the Laplacian where w q is used as weighting factors that optimise the accuracy of the approximation [7].

B. Boundary conditions of the lattice-Boltzmann algorithm
The boundary conditions for the numerical method are such to give the correct representation of the diffuse interface system.For the solid boundary conditions, we implement an interpolated bounceback algorithm [8], for the unknown particle population f q travelling in the opposite direction to the solid surface, that is, q is such that c q +c q = 0.The distance fraction, δ q , is such that where x b is the position of the solid surface.In Equation ( 8) f ⋆ q corresponds to the prestreaming particle population, i.e., before updating the populations to t + 1.This ensures that the velocity of the boundary is as in Eq. ( 12) of the manuscript and conserves mass density as in Eq. ( 13) of the main text.Similarly, ensures that the phase-field is conserved and no diffusive fluxes permeate the solid bound-

aries.
At the open boundary, we implement an anti-bounceback algorithm [9].Then, which satisfies Equation (17) of the main text and Equation ( 16) by extrapolating the velocity at the boundary, u i = u(x), which is assumed to be at x + c q (1/2).Similarly, which imposes the value of the phase field and chemical potential at the boundaries as in Eq. ( 14) and (15) of the manuscript, respectively.
For further details and validation of the numerical method and the different boundary conditions see Ref. [10].

II. ADDITIONAL SIMULATION RESULTS
A. Chemical potential gradients

Figure 2 FIG. 2 :Figure 3
Figure2shows the surface tension profile of the strongest condensation gradient we explored.The change in surface tension along the interface of the droplet is negligible up to the point where the droplet approaches the boundary, generating stronger gradients due to the imbalance caused by the boundary conditions.Before this point, the apparent motion of the droplet is predominantly dominated by the phase change, as seen in Figure1d).